#include <fstream>
#include <iostream>
#include <iomanip>
#include <stdio.h>
#include "TFile.h"
#include "TTree.h"
#include "TLine.h"
#include "TMath.h"
#include "TSystem.h"
#include "TStyle.h"
#include "TGraph.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TEllipse.h"
#include "TClonesArray.h"

using namespace std;

float acc[5][305*34]={{0}};
int    nphi[15] = {3,5,8,11,14,17,20,24, 27, 30, 33, 36, 39, 19, 19}; //305
int sumnphi[15] = {0,3,8,16,27,41,58,78,102,129,159,192,228,267,286}; //305


double Acc(int pmtno,float xpos,float ypos,float zpos);
void test(){

  int ret = gSystem->Load( "../Objects/libObjectClass.so" );
  if( ret!=0 ){
    cout << "Failed to Load Library: ret = " << ret << endl;
    return;
  }
  
  TFile *ff=new TFile( "/scratch/shinichi/OutputFiles/75L_RPhi/75L_Full_2phase.root" );
  TTree *tree=(TTree*)ff->Get("tree");

  int nPhoton,nPMThits[5];
  tree->SetBranchAddress("nPhoton",&nPhoton);
  tree->SetBranchAddress("nJPMThits",&nPMThits[0]);
  tree->SetBranchAddress("nMPMT1hits",&nPMThits[1]);
  tree->SetBranchAddress("nMPMT2hits",&nPMThits[2]);
  tree->SetBranchAddress("nMPMT3hits",&nPMThits[3]);
  tree->SetBranchAddress("nMPMT4hits",&nPMThits[4]);

  int nev = tree->GetEntries();
  cout<<"NEV = "<<nev<<endl;

  int NPHOTON=0;
  int NPMTHITS[5]={0};
  int IEV = 0;

  for(int iev=0;iev<nev;iev++){

    tree->GetEntry(iev);

    NPHOTON += nPhoton;
    for(int ipmt=0;ipmt<5;ipmt++) NPMTHITS[ipmt] += nPMThits[ipmt];
    
    if((iev%10)==9){

      int runid = IEV;
      int id=-1;
      int runID= runid % 305;
      if( runID>=0  && runID<sumnphi[1] )  id=0;
      else if( runID<sumnphi[2] )  id=1;
      else if( runID<sumnphi[3] )  id=2;
      else if( runID<sumnphi[4] )  id=3;
      else if( runID<sumnphi[5] )  id=4;
      else if( runID<sumnphi[6] )  id=5;
      else if( runID<sumnphi[7] )  id=6;
      else if( runID<sumnphi[8] )  id=7;
      else if( runID<sumnphi[9] )  id=8;
      else if( runID<sumnphi[10])  id=9;
      else if( runID<sumnphi[11])  id=10;
      else if( runID<sumnphi[12])  id=11;
      else if( runID<sumnphi[13])  id=12;
      else if( runID<sumnphi[14])  id=13;
      else if( runID<305 )         id=14;
      else if( id==-1 ){ cout<<runid<<endl;return;}
      
      int iphi = (runID-sumnphi[id]) % nphi[id];
      int ir   = id;
      int iz   = runid / 305 ;
      
      float rr = 5.+10.*ir;
      float phi = iphi*(TMath::Pi()/(float)(nphi[ir]-1)); //rad
      float zz;
      if( iz<24 ) zz = -595.+iz*10.;
      else        zz = -360.+iz*1.;

      for(int ipmt=0;ipmt<5;ipmt++){
	acc[ipmt][IEV] = NPMTHITS[ipmt]/(float)NPHOTON;
      }

      NPHOTON = 0;
      for(int ipmt=0;ipmt<5;ipmt++) NPMTHITS[ipmt] = 0;
      IEV++;
    }
  }

  float x=32;
  float y=0;
  float z=-465;

  TH2F *test[5];
  TH2F *test2[5];
  TH1F *test3[5];
  TGraph *g1=new TGraph(200);

  for(int i=0;i<5;i++){  
    test[i]=new TH2F(Form("htest_%d",i),Form("htest_%d",i),300,-150,150,300,-150,150);
    test2[i]=new TH2F(Form("htest2_%d",i),Form("htest2_%d",i),300,-150,150,300,-620,-320);
    test3[i]=new TH1F(Form("htest3_%d",i),Form("htest3_%d",i),3000,-150,150);
  }
  // for(int ix=-150;ix<150;ix++){
  //   for(int iy=-150;iy<150;iy++){
  //     for(int i=0;i<5;i++){
  // 	test[i]->Fill(ix,iy,Acc(i,ix,iy,z));
  //     }
  //   }
  // }

  // for(int i=0;i<5;i++){
  //   for(int ix=-150;ix<150;ix++){
  //     for(int iz=-600;iz<-350;iz++){
  // 	test2[i]->Fill(ix,iz,Acc(i,ix,y,iz));
  //     }
  //   }
  // }
  for(int i=0;i<5;i++){
    for(int ix=-100;ix<100;ix++){
      test3[i]->Fill(ix/10.,Acc(i,ix/10.,y,z));
      if(i==0){
	g1->SetPoint(ix+100,ix/10.,Acc(i,ix/10.,y,z));
	cout<<ix/10.<<" "<<Acc(i,ix/10.,y,z)<<endl;
      }
    }
  }

  float max = 0.06;

  gStyle->SetOptStat(000000);
  int pos[5]={5,3,1,7,9};
  char cpmt[5][64]={"JPMT","Meg5","Meg3","Meg4","Meg1"};

  TCanvas *c1=new TCanvas("c1","c1",1200,1000);
  c1->Divide(3,3);
  for(int i=0;i<5;i++){
    c1->cd(pos[i]);
    char name[128];
    sprintf(name,"%s | Z = %d (mm)",cpmt[i],(int)z);
    test[i]->SetTitle(name);
    test[i]->SetXTitle("X (mm)");
    test[i]->SetYTitle("Y (mm)");
    test[i]->SetMaximum(max);
    test[i]->Draw("colz");
  }

  TCanvas *c2=new TCanvas("c2","c2",1200,1000);
  c2->Divide(3,3);
  for(int i=0;i<5;i++){
    c2->cd(pos[i]);
    char name[128];
    sprintf(name,"%s | Y = %d (mm)",cpmt[i],(int)y);
    test2[i]->SetTitle(name);
    test2[i]->SetXTitle("X (mm)");
    test2[i]->SetYTitle("Z (mm)");
    test2[i]->SetMaximum(max);
    test2[i]->Draw("colz");
  }

  TCanvas *c3=new TCanvas("c3","c3",1200,1000);
  c3->Divide(3,3);
  for(int i=0;i<5;i++){
    c3->cd(pos[i]);
    char name[128];
    sprintf(name,"%s | Y = %d , Z = %d (mm)",cpmt[i],(int)y,(int)z);
    test3[i]->SetTitle(name);
    test3[i]->SetXTitle("X (mm)");
    test3[i]->SetYTitle("Acceptance");
    test3[i]->SetMaximum(max);
    test3[i]->Draw("hist");
  }


  TCanvas *c4=new TCanvas("c4","c4",1200,1000);
  g1->Draw("AP");


  // c1->SaveAs(Form("./plot20130309/2phase_hXY_z%d.png",(int)z));
  // c2->SaveAs(Form("./plot20130309/2phase_hXZ_y%d.png",(int)y));

}

//----- Function
// -150 < xpos,ypos < 150 (mm)
// -600 < zpos < -350;
// pmtno=0; J-PMT
// pmtno=1; Meg5
// pmtno=2; Meg3
// pmtno=3; Meg4
// pmtno=4; Meg1

double Acc(int pmtno,float xpos,float ypos,float zpos){

  if(pmtno>4){cout<<"--- Wrong Input PMT No. ---"<<endl;return 0;};
  if( zpos>=-350.||zpos<=-600. ) return 0;
  
  float RR,theta;
  int ir,ir2;
  
  RR   = sqrt(xpos*xpos+ypos*ypos);
  if( RR>=150.) return 0.;
  
  if( RR<=0. ) RR=0.01;    
  theta = TMath::ACos( xpos/RR );
  
  //-----
  ir = (int)(RR/10.);
  if(ir>14) return 0;

  float dphi = TMath::Pi()/(float)(nphi[ir]-1);
  float dphi2;
  if( ir==(int)(RR/10.-0.5) ){
    if( ir==14 )  ir2 = ir;
    else          ir2 = ir+1;
    dphi2 = TMath::Pi()/(float)(nphi[ir2]-1);
  }
  else{
    if( ir==0 )	ir2 = ir;
    else          ir2 = ir-1;
    dphi2 = TMath::Pi()/(float)(nphi[ir2]-1);
  }
  //-----
  int no1=-1;
  for(int j=0;j<nphi[ir];j++){
    float tmpphi = j*dphi;
    float delta = theta-tmpphi;
    if( delta<0 ){ no1=j; break;}
  }
  if( no1<0 ) no1=nphi[ir]-1;
  
  int no2=-1;
  for(int j=0;j<nphi[ir2];j++){
    float tmpphi = j*dphi2;
    float delta = theta-tmpphi;
    if( delta<0 ){ no2=j; break;}
  }
  if( no2<0 ) no2=nphi[ir2]-1;
  
  //-----
  int iza,izb;
  float dz,zza,zzb;
  if(zpos<-360){
    iza = (int)((zpos+600)/10);
    if(iza*10-595<zpos){
      izb = iza+1;
      if(iza==23) izb=iza;
    }
    else if(zpos<=-595)   izb = iza;
    else                  izb = iza-1;  
    dz = 10.*(iza-izb);
    zza = -595.+iza*10.;
    zzb = -595.+izb*10.;
  }
  else{
    iza = (int)(24+zpos+360);
    if(iza==33) izb = iza;
    else        izb = iza+1;
    dz = 1.;
    zza = -360.+(iza-24)*1.;
    zzb = -360.+(izb-24)*1.;
  }
  
  //-------
  if(ypos<0&&pmtno!=0) pmtno = 5-pmtno;
    
  float Acc1_1a = acc[pmtno][305*iza+sumnphi[ir]+no1];
  float Acc1_1b = acc[pmtno][305*izb+sumnphi[ir]+no1];
  float Acc1_2a = acc[pmtno][305*iza+sumnphi[ir]+no1-1];
  float Acc1_2b = acc[pmtno][305*izb+sumnphi[ir]+no1-1];
  float Acc2_1a = acc[pmtno][305*iza+sumnphi[ir2]+no2];
  float Acc2_1b = acc[pmtno][305*izb+sumnphi[ir2]+no2];
  float Acc2_2a = acc[pmtno][305*iza+sumnphi[ir2]+no2-1];
  float Acc2_2b = acc[pmtno][305*izb+sumnphi[ir2]+no2-1];
  
  float Acc1a = (Acc1_1a-Acc1_2a)/ dphi*theta + (no1* dphi*Acc1_2a-(no1-1)* dphi*Acc1_1a)/dphi ;
  float Acc1b = (Acc1_1b-Acc1_2b)/ dphi*theta + (no1* dphi*Acc1_2b-(no1-1)* dphi*Acc1_1b)/dphi ;
  float Acc2a = (Acc2_1a-Acc2_2a)/dphi2*theta + (no2*dphi2*Acc2_2a-(no2-1)*dphi2*Acc2_1a)/dphi2 ;
  float Acc2b = (Acc2_1b-Acc2_2b)/dphi2*theta + (no2*dphi2*Acc2_2b-(no2-1)*dphi2*Acc2_1b)/dphi2 ;
  if(Acc1a<0) Acc1a=0;
  if(Acc1b<0) Acc1b=0;
  if(Acc2a<0) Acc2a=0;
  if(Acc2b<0) Acc2b=0;
  
  float Acca,Accb;	
  if( ir==ir2 ){
    Acca = Acc1a; 
    Accb = Acc1b; 
  }
  else{
    Acca = (Acc1a-Acc2a)/(10.*(ir-ir2))*RR + ((10.*ir+5.)*Acc2a-(10.*ir2+5.)*Acc1a)/(10.*(ir-ir2)) ;
    Accb = (Acc1b-Acc2b)/(10.*(ir-ir2))*RR + ((10.*ir+5.)*Acc2b-(10.*ir2+5.)*Acc1b)/(10.*(ir-ir2)) ;
  }
  if(Acca<0) Acca=0;
  if(Accb<0) Accb=0;

  double Acc;
  if(zza==zzb)    Acc = Acca;
  else            Acc = (Acca-Accb)/(zza-zzb)*zpos + (zza*Accb-zzb*Acca)/(zza-zzb);
  if(Acc<0) Acc=0;

  if(pmtno==0&&zpos==-465.&&ypos==0.&&fabs(xpos)<=10.) cout<<"xpos="<<xpos<<" "<<ir<<" "<<ir2<<" "<<no1<<" "<<no2<<" "<<Acc<<endl;
  
 
  return Acc;

}

